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TRAPEZOIDAL  MONTE  CARLO  INTEGRATION* 

ELIAS  MASRYt  and  STAMATIS  CAMBANISt 

Abstract.  The  approximation  of  weighted  integrals  of  random  processes  by  the  trapezoidal  rule  based 
on  an  ordered  random  sample  is  considered.  For  processes  that  are  once  mean-square  continuously 
differentiable  and  for  weight  functions  that  are  twice  continuously  differentiable,  it  is  shown  that  the  rate 
of  convergence  of  the  mean-square  integral  approximation  error  is  precisely  n~4,  and  the  asymptotic  constant 
is  also  determined. 

Key  words.  Monte  Carlo  integration  of  random  processes,  trapezoidal  rule,  rate  of  quadratic-mean 
convergence 

AMS(MOS)  subject  classifications.  65D30,  60G1 2,  65U05,  62G0S 

1.  Introduction,  results,  and  discussion.  The  simplest  Monte  Carlo  method  for 
approximating  the  integral 

0.1)  /(/)=[/(»)* 

Jo 

of  a  (square  integrable)  function  /  over  a  finite  interval,  uses  n  independent  samples 
t/| ,  •  •  • ,  U„  from  a  uniform  distribution  over  the  unit  interval  and  forms  the  average 
estimate 

(1.2)  =  -  i  f{U,). 

n  (_ i 

J„  is  an  unbiased  estimator  of  /:  EJn(f)  ~  1(f),  and,  in  view  of  the  strong  law  of  large 
numbers,  it  is  consistent,  i.e.,  as  the  sample  size  n  increases  to  infinity,  J„(f)  tends  to 
1(f)  with  probability  one,  i.e.,  for  almost  every  realization  of  The  variance 

or  mean-square  error  of  J„(f)  is 

(1.3)  E  [1(f)  -Jn(f)  J2  = 1  {1(f)  -V(f)Y) 

n 

and  thus  tends  to  zero  at  the  rate  of  n~'.  As  the  constant  1(f)  -  l7(f)  is  strictly  positive 
for  all  functions  /  which  are  not  almost  everywhere  constant,  no  improvement  in  the 
rate  of  convergence  is  to  be  expected  from  any  smoothness  assumptions  on  f 

Yakowitz,  Krimmel,  and  Szidarovszky  [9]  proposed  improving  the  convergence 
rate  in  (1.3)  of  the  crude  Monte  Carlo  method  by  using  quadrature  formulas  instead 
of  the  simple  averaging  in  (1.2).  They  specifically  studied  the  trapezoidal  rule,  based 
on  the  ordered  sample  o^0<  fn>,  <  tn_2<  •  •  •  <  1  A  fn(1+1  obtained  from  the  indepen¬ 
dent,  uniformly  distributed  samples  l/,,  U2,  •  •  • ,  U„,  which  is  given  by 

(1.4a)  Uf)  =  \  1  [/(L.<)+/(L,+,)](L.i+,-L,) 
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and  uses  n  4-2  sample  points,  the  n  random  samples,  and  the  fixed-interval  endpoints. 
Since  the  trapezoidal  rule  can  be  written  in  the  form 

(l.4b)  /„(/)=^{^,/(0)+  i  (t „,l+1  - +(l  - r„.J/(l)} 

it  can  be  considered  as  a  weighted  Monte  Carlo  rule  with  random  weights.  When  / 
has  continuous  second  derivative  they  proved  that 

(1-5)  E[/(f)-I„(f)]2S^fi 

n 

for  some  constant  C(f )  and  all  nil,  and  they  provided  simulation  evidence  that  the 
convergence  rate  of  the  mean-square  error  is  in  fact  n~ \ 

In  this  paper,  we  consider  weighted  integrals  of  random  processes  and  establish 
the  rate  of  convergence  and  asymptotic  constant  for  the  trapezoidal  rule  (1.4a)  based 
on  an  ordered  random  sample.  As  a  consequence,  the  rate  of  convergence  and 
asymptotic  constant  for  integrals  of  nonrandom  functions  are  also  determined. 

Throughout  X  =  {X(t,  w),  0^  t  S  1}  is  a  second-order  random  process  with  mean 
zero,  EX(t)  =  0,  and  covariance  function  R(t,s)  =  E{X(t)X(s)},  defined  on  some 
probability  space.  We  assume,  with  no  further  notice,  that 

(i)  The  random  process  X  and  the  random  sample  { L/,}  are  mutually  independent, 

and 

(ii)  X(t,a>)  is  jointly  measurable  in  t  and  to. 

The  first  assumption  simplifies  the  task  of  computing  expected  values  and  is  quite 
natural  as  the  randomness  of  the  sampling  mechanism  is  generally  in  no  way  dependent 
on  or  related  to  the  randomness  of  the  stochastic  integrand.  The  second  assumption 
enables  us  to  integrate  with  respect  to  t  for  almost  every  fixed  path  to  and  to  obtain 
as  a  result  of  the  integration  a  random  variable,  whose  expectation  may  therefore  be 
computed.  This  is  a  .minimal  regularity  assumption  on  the  process  X;  when  X  is 
mean-square  continuous,  as  we  will  shortly  assume,  it  always  has  a  jointly  measurable 
version.  Following  standard  practice,  we  delete  the  probability  variable  to  and  write 
X(t)forX(r,  to)  (just  as  we  write  tni  and  t/,  for  the  random  variables  tni(to)  and  Ufto)). 

We  first  derive  under  general  conditions  an  explicit  expression  for  the  mean-square 
error  in  the  approximation  of  the  random  integral: 

/(/*)=  \'f(t)X(t)dt 
Jo 

by  the  trapezoidal  rule  based  on  ordered  random  samples: 

In(fX)=\  I  [/(fB,i)X(tB,()+/(tn,i+I)X(fn,+  I)](fn,(+J-fn,(). 

I  (-0 

This  expression  can  be  used  to  evaluate  finite  sample  size  performance. 

Theorem  1.  Ifi'0R(t,  t)f\t)  dt  <  oo,  then  for  all  n  e  1  we  have 

E[I(fX)  -  I„(fX)]2 

(1.6a)  fi(0,0)/2(0)  +  *(l,  l)/2(l)  +  fi(0,  l)/(0)/(l) 

2(n  +  l)(n  +  2) 
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(1.6b) 


['{«(/,  mi)f(t)  +  R(0, 1  -  t)/(0)/(l- 1)} 

Jo 

f  (l-0n+I  n2  +  3n -2  n+l  n  +  1  „  n  ...l  . 

1-Wl7+^TTT'  —  ,  +i'  r' 

(,«c 

ff  R(t,  s)f(t)f(s) 

J  J 05f <s<l 


3  + --  [t"+l  +  (1  - 0"+l]-~ [/"  +  (1  - *)"]]  dt 


(1.6d)  + 


[r  +  (l-Sn-i(s-/)"+i«(n-l)(l-s  +  /)"-2 

2  4 


n2(1-s+l)""'+J  (H-2HH  +  3HI  -»+ 


'I 


Using  the  expression  of  the  mean-square  error  in  (1.6),  we  now  show  that  when 
the  function  /has  two  continuous  derivatives  and  the  process  X  has  one  quadratic-mean 
derivative  which  is  mean-square  continuous,  then  the  rate  of  convergence  of  the 
mean-square  integral  approximation  error  is  precisely  n~4;  and  we  also  determine  the 
asymptotic  constant. 

Theorem  2.  If  f(t)  has  continuous  second  derivative  on  [0,1]  and  R(t,s)  has 
continuous  mixed  partial  derivatives  R‘~J(t,  s)  of  order  2, 0  S  /  +j  S  2,  on  the  unit  square 
[0, 1]  x  [0, 1]  and  of  order  3,  i+j  =  3,  off  its  diagonal,  then 


(1-7) 


£[/(/*) -/„(/X)]2  = 


C{f  J?)  +  o(l) 

(n  +  l)(n  +  2)(n  +  3)(n  +4)' 


The  asymptotic  constant  in  (1.7)  is  given  by 
2C(/  K)=i/?(0,0)/(0)2  +  iJ?(l,  l)f(0)f(l) 

+  /?■>,  0)/(0)/(0)  +  K,  o(l,  l)/(l)/(l)-/?’  °(0,  l)/(0)/'(l) 

(1.8a)  - RO-l(0,  l)/(0)/(l)  +  3J?2-°(l,  1)/2(1) -3/?2O(0, 0)/2(0) 

+  2R (0, 0)/2(0)  -  R '■'(  1, 1  )/2(  1 )  -  R (0, 1 )/( 0)/(  1 ) 

+  3  f  [/?***( f,  t)-2R2’°(t,  t)]f{t)f(t)  dt-3  R3fi(t-,t)f\t)dt. 


By  putting  in  Theorem  2,  R(t,  s)  =  l,  we  obtain  the  r*rp  ise  rate  conjectured  in 
[9]  for  the  integral  approximation  of  nonrandom  funu.on  ,  plus  of  course  the 
asymptotic  constant. 

Corollary.  If  f  has  continuous  second  derivative  on  [0, 1]  then 


(1.9a) 

Thus, 

(1.9b) 


E[I(f)-I„(f)T  = 


[fw-rm2+o{i) 

4(n  +  l)(n  +  2)(n  +  3)(n  +  4)‘ 


lim  n4E[I(f)  -  I„(f)Y  =  it/d) ~/(0)]2. 

n-*oo 


It  is  interesting  to  note  that  the  asymptotic  constant  in  (1.9b)  for  the  trapezoidal  rule 
with  random  samples  has  the  same  functional  form  as  the  classical  asymptotic  constant 
for  the  trapezoidal  rule  with  equally-spaced  samples  [7,  Thm.  3.3],  but  is  larger  by  a 
factor  of  36. 
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Returning  to  the  general  case  of  Theorem  2,  we  note  that  the  assumption  of 
continuous  mixed  partial  derivatives  of  R  of  order  up  to  2  on  (0. 1  ]:  is  equivalent  to 
the  assumption  that  X  has  one  mean-square  continuous  quadratic-mean  derivative. 
The  additional  assumption  of  differentiability  of  order  3  of!  the  diagonal  is  weak,  and 
is  always  satisfied,  e.g.,  when  X  is  stationary,  has  rational  spectral  density,  and  exactly 
one  quadratic-mean  derivative. 

The  expression  of  C(/  R)  in  (1.8a)  is  not  symmetric  and  cannot  be  symmetrized 
under  the  conditions  of  Theorem  2.  If  R(i,s)  has  continuous  mixed  partial  derivatives 
of  order  3  throughout  the  unit  square  (rather  than  merely  of!  its  diagonal)  then  the 
asymptotic  constant  in  (1.8a)  takes  the  following  simple  and  symmetric  form: 

(1.10)  C,ym(/,R)  =  ‘{M"(0,0)  +  M,-,(1.1)-2M"(0.1)|. 


where  M(t,  s)  =/(/)/?(/,  s)f(s).  The  expression  within  braces  in  (1.10)  equals 
£[(/*)'( l)-(/*)'(0)]J,  where  prime  denotes  quadratic-mean  derivative.  Thus  under 
the  slightly  more  stringent  assumption  of  continuous  mixed  partial  derivatives  of  order 
3  on  the  unit  square  we  obtain 


(1.11) 


E[l(fX)-IH(fX)]2  = 


EUfxn\)-(fxno))2+o(\) 

4(n  +  l)(n  +  2)(«  +  3)(n  +  4)  ' 


This  expression  is  the  stochastic  analogue  of  the  nonrandom  case  in  (1.9a)  and  it  holds 
even  though  X  is  not  assumed  to  have  a  second  quadratic-mean  derivative.  It  is  dear 
from  (1.11)  that  in  general  the  asymptotic  constant  is  positive  and  no  faster  rate  of 
convergence  can  be  achieved  by  requiring  any  further  smoothness  of  /  or  of  R.  Under 
the  conditions  of  Theorem  2,  the  asymptotic  constant  differs  from  its  symmetric 
expression  by  the  following: 


C(/,K)  =  Ciym(/,/?)+| 
(1.8b) 


[*'•'(/,  f)-2K3-°(f,r)]/(r)/’(r)df 


-  j  («'•'('.  ‘)-2R20(t,t)f2(t)]l 


R’°(»-,r )/-’(/)  dt 


) 


When  the  random  process  X  is  stationary,  i.e.,  R(t,  s)  =  R(t  -  s ),  then  the  general 
form  (1.8a)  of  the  asymptotic  constant  simplifies  to 


2 CJf,  R)  = X-  R(0)Lr(0)2  +/( 1  )2]  -  R(  1  1 )  +  R’{  1  )[/(0)/(  1 )  -/( 0)/(  1 )] 


(1.12)  R"(0)[/2(0)+/2(l)]  +  /?"(l)/(0)/(l)-3R"(0-)  J*  f2 

=  l-  E[(fX)'(  1)  -  (PO'(O)]2  -  3iT(0-)  | '  f2. 

When  instead  of  using  in  the  trapezoidal  rule  (1.4a)  ordered  random  samples,  we 
use  equidistant  samples  f,  =  i/(n  + 1),  i  =  0,  1,  •  •  • ,  n  +  1,  then  it  has  been  shown  in 
[2,  App.  B]  that  for  stationary  processes,  under  the  assumptions  of  Theorem  2,  we  have 

E[l(fX)  -  I„{fX)]2  =  CgtA*)*0*1) 
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where 


CM  R )  =  ^ R mno)2 +  a  1  )2] -  R{  1  )/(0)/(  1 )  +  R '( 1  )[/(0)/(  1 )  -/ (0)/(  1 )] 

-  \  RH(0)[f2(0)  +/2( \)]+R"(l  )/(0)/(  1 )}  -  “  R"'(  0-)  J 1  f2 

or 

CM R)=t^£[(/a:)'(i)-(/x)'(o)]2+-^[-r"'(o-)]  f/2<^ cm, R). 

144  360  Jo  36 

It  follows  that  (CJ Ced)l/4>  (36) 1/4  ^  2.45  and  thus,  asymptotically,  at  least  two-and-a- 
half  times  more  random  samples  are  required  than  equidistant  samples  for  the  same 
accuracy  measured  in  mean-square  error.  It  is  of  course  quite  natural  that  equidistant 
samples  provide  a  superior  approximation  than  ordered  random  samples  with  the  same 
average  distance. 

The  analysis  carried  out  here  suggests  that  kth-order  quadrature  rules  based  on 
ordered  random  samples  should  have  mean-square  error  with  rate  of  convergence 
„-«*+■>  when  acting  on  nonrandom  functions  with  continuous  (fc  +  l)st  derivative,  or 
random  processes  with  (essentially)  mean-square  continuous  kth  quadratic-mean 
derivative,  or  on  their  products.  (The  trapezoidal  rule  considered  here  and  in  [9]  is  a 
first-order  quadrature  rule.) 

It  should  be  finally  mentioned  that,  for  integrals  of  nonrandom  functions,  Haber 
has  developed  a  stratified  Monte  Carlo  rule  with  rate  n~ 3  [3];  a  stratified  and  sym¬ 
metrized  Monte  Carlo  rule  with  rate  n~s  [4];  and  certain  stratified  stochastic  quadrature 
formulas  with  rate  n~'~2k  when  approximating  the  integral  of  a  nonrandom  function 
/  with  continuous  fcth  derivative  [5].  For  weighted  integrals  of  random  processes,  a 
simple  Monte  Carlo  rule  with  rate  n  '  and  a  stratified  Monte  Carlo  rule  with  rate  n~3 
have  been  developed  by  Schoenfelder  [6]  (see  also  [1]). 

Example.  We  illustrate  by  an  example  the  finite  sample  size  performance  of  the 
trapezoidal  rule  with  random  sampling.  We  consider  the  stationary  process  X(t)  with 
correlation  function 


R(t)  =  (1  +  )3|f|)  e~3|r| 

where  p  >  0.  Note  that  R(0)  =  1  and  that  X  has  precisely  one  quadratic-mean  deriva¬ 
tive.  For  simplicity  we  take  /(r)“l  so  that  the  random  integral  to  be  estimated  is 
/(X)  =  Jo  X(t)  dt  and  its  estimate  is  /„(X)  of  (1.4a).  The  variance  a2  of  /(X)  is  given 
by 


Using  Theorem  1,  we  find  after  some  algebra  that  the  mean-square  error  is  given  by 

-3  +  4n2  +  n  +  2) 


£[/(X)-Jb(X)] 


1  ^(n  +  3)(n: 


2/82 


rt4+4n3+6n2  +  9n  +  2  n*+4n3+5n2-4n-4) 


0(n  + 1) 


-+- 


2(n  +  l)(fi  +  2) 
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(1-13) 


for  nS  1,  where 


i  -/if  8  |  2w  +  5  t  («  +  3): 

+  2e  l(n  +  l)(n  +  2)+(n  +  l)(n  +  2)  /3(n  +  l) 

+  ^a(fi,i8)|(n  +  2)-^(n  +  3)-} 

+~o(n  -  1,  -f})n  e~p  j/}(n  +  l)  +  (3/r +  6n  +  5) 


+  —  (3/iJ  +  12n2+  13/i-HO) 

8 


4 — 5  (**  +  3)(«*'+4/r  + 


n  +  2)J 


a(n,/3)= 

Jo 


x"dx. 


For  n  =0,  the  mean-square  error  can  be  computed  directly  yielding 


,1  6  R  [ 5  8  6  6 

E[HX)-I0(X)f=rj2+e-'[-+Z+-+-2\. 


The  asymptotic  constant  C„(  1,  R)  is  given  by 


C,t(l,R)^{l  +  68  +  (8-\)e^}~r 


Let  m  -  2,  3,  ■  •  •  be  the  (true)  sample  size,  m  =  n  +  2,  with  corresponding  mean- 
square  error  mse  ( m )  =  E[I(X)  -  /m_2(X)]2.  The  fractional  mean-square  error  is  then 
given  by  mse  (m)/<j2. 

In  selecting  appropriate  values  of  8  for  numerical  display  of  the  finite  sample 
size  performance,  the  behavior  of  the  fractional  error  mse  (0)/ a-2  (based  only  on  the 
endpoints  X(0)  and  X(l))  as  a  function  of  /3  was  investigated.  It  is  seen  from  Table 
1  that  for  values  8  =  1  the  fractional  mean-square  error  is  too  small  so  that  10{X) 
already  provides  a  fairly  accurate  approximation  of  I{X).  We  select  therefore  two 


Table  1 

The  fractional  error  mse  (0 )/cr2  and  the  asymptotic 
constant  C„(1,  R)  as  functions  of  f). 


p 

mse  (0)/ir2 

C„(1,K) 

0.2 

2.3624  x  10'4 

3.09  xIO'2 

0.4 

1.686  x  10-3 

0.24 

1 

1.929  xlO-2 

3.5 

2 

9.862  x  I0~2 

26.27 

3 

0.225 

85.95 

4 

0.377 

200.44 

5 

0.5376 

387.84 

6 

0.698 

666.22 

7 

0.854 

1053.6 

8 

1.0058 

1568.08 

9 

1.152 

2227.54 

10 

1.295 

3050.02 
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values  p  =  3  and  p  =  5  corresponding  to  moderate  values  of  mse  (0)/<r2.  The  asymptotic 
constant  Csl(l,  R)  is  monotonically  increasing  with  p  and  is  asymptotically  equal  to 
3/33  for  large  p.  Table  1  provides  a  few  typical  values. 

In  Fig.  1  the  fractional  mean-square  error  mse  (m)/<r2  is  plotted  as  a  function  of 
the  sample  size  m  =  2,  •  ■  • ,  26,  for  P  =  3  and  p  =  5.  It  is  seen  that  for  the  smaller  value 
of  p  =  3,  the  fractional  error  is  considerably  smaller  for  each  sample  size  m.  This  can 
be  explained  by  the  less  rapid  decay  of  R(t)  and  hence  the  larger  correlation,  on  the 
average,  between  consecutive  samples  so  that  I„(X)  provides  a  better  estimate  of  I(X) 
in  this  case.  The  closeness  of  the  fractional  mean-square  error  to  its  asymptotic  value, 


mse  (m)/<r2~ 


CJX  R)/ <t2 
(m  -  \  )m(m  +  l)(m  +  2)’ 


is  displayed  in  Fig.  2  for  parameter  p  =  3.  Note  that  the  asymptotic  value  overestimates 
the  true  error  for  all  sample  sizes  m  in  the  plotted  range.  Naturally  the  discrepancy 
between  the  two  values  diminishes  as  m  increases.  To  see  more  clearly  the  convergence 
of  the  scaled  mean-square  error  (m  -\)m(m  +  l)(m  +  2)  mse  (m)  to  the  asymptotic 
constant  C„(l,  R)  as  m  increases,  we  display  in  Fig.  3  the  values  of  these  two  quantities 
for  m  =  2,  •  •  • ,  26  with  the  chosen  values  of  p  as  parameter.  Again  for  the  smaller 
p  -  3  the  discrepancy  between  these  two  quantities  is  smaller  for  each  m  =  2,  •  •  • ,  26 
than  for  >3  =  5. 

2.  The  mean-square  error.  In  this  section  we  derive  the  exact  expression  of  the 
mean-square  error  given  in  Theorem  1. 

Since  the  trapezoidal  rule  approximation  /„  to  the  integral  I  is  based  on  the 
ordered  samples  tn  X  <  t„  2  <  •  •  •  <  t„  „  obtained  from  n  independent  uniformly  dis¬ 
tributed  samples  on  (0, 1 ),  we  need  certain  properties  of  the  order  statistics  from  the 
uniform  distribution,  which  we  summarize  first.  For  brevity  we  will  write  tk  for  tn  k. 


Fig.  1.  Fractional  mean-square  error  mse  (m)/o2  as  a  function  of  the  sample  size  m. 


30 


Fractional  Mean-Square  Error 


232 


ELiAS  MASRY  AND  STAMATIS  C  AMBANIS 


0  10  20  30 

Sample  Size 


Fig.  2.  Exact  and  asymptotic  fractional  mean-square  error  as  functions  of  the  sample  size  (/3  =  3). 


Fig.  3.  Scaled  mean-square  error  (m  -  l)m(m  +  l)(m  +  2)  mse  (m)  as  a  function  of  the  sample  size  m. 
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The  joint  distribution  of  order  statistics  is  an  ordered  multivariate  Dirichlet 
distribution;  see  Wilks  [8,  §§  8.7.5  and  8.7.7].  Specifically,  the  ordered  samples  u, . 
tk,+k,,  -M+k2+  -+k„,  (where  the  k,' s  are  positive  integers  with  k,  +  k2  +  •  ■  -  +  kmSn)  have 
joint  probability  density  function  denoted  by  />*,.*,+*,...*,+• +kjxu  x2,  ■  •  ■  ,  xm)  and 
given  by 

_ r(n  +  l) _ 

r(k,)-r(fcm)r(n+i-/c1 - km > 

•  •  •  •  (xm-xm_l)‘»-'(l -Jt.)""*'""”4- 

for  0<x, <x2<-  •  -<xm<l  and  zero  elsewhere.  We  will  make  explicit  use  of  the 
following  expressions. 

The  smallest-  and  largest-order  statistics  /,  and  ln  have  densities 

(2.1)  Pi(x)  =  n(l  -x)"-1,  pn(x)  =  nx"~',  0<x<l, 

and  joint  density 

(2.2)  p,.n(x,y)-n(n-l)(y-x)n~2,  0<x<y<l. 

Two  consecutive  order  statistics  f,  and  t,+,  (1  S  i  S  n  -  1)  have  joint  density 

xi_l  (i —  v)n_l_‘ 

(2.3)  Pu+i(x,y)  =  nl— —  (n_1_i), ,  0<x<y<l. 

We  will  need  the  value  of  their  sum 

»- 1  n—2  (1  _  y)»-2  * 

(2.4)  X  pu+l(x,y)  =  n\  X  —  -  -  -  ~  ~  =  n(n  -  1)(1  -y  +  x)"~;. 

/=i  k-o  a!  (n  —  z  —  k)l 

Three  consecutive  order  statistics  th  ti+l,  t,+2  (1  SiSn-2)  have  joint  density 

x-'  (1  -z)n_2“' 

Pu+u+zU y,*Y =  » !  (7T7)j  („_2 _ •  0<x<y<z<  1, 

and  their  sum  will  be  used: 

n-2  n-2  X*  (1  ~Z)"~3~k 

(2.5)  X  Pi,i+\.i*2(x, y,  z)  =  n\  X  7:7 - - r77  =  n(n  -  l)(n  —  2)(  1  -  z  +  x)"  \ 

/= 1  k  =0  a'  (n  -3  -  fcj! 

We  will  also  use  the  following  two  trivariate  densities 

.  v  (y-x)1-2  (1  -z)"-1-1 

Pi.u+i(x,y,z)  =  n\  —7 - 7—777,  0<x<y<z<l,  2SiSn-\, 

(i-2)!  (n-l-ij! 


,  .  x'-'  (z-y)"-2- 


0<x<y  <  z  <  1,  IS  iS  n-2. 


and  their  sums 


I  Pi.i,H-iU3'»z)  =  n(w“1)("-2)(l-z+y-x)"  \ 

i”  2 


I  Pi,i+\.n(x,  y,  z)  =  n(n  -  l)(n-2)(z-y  +  x)n  \ 


Finally  two  pairs  of  consecutive  order  statistics  rf,  tJt  tJ+,  (l<i+l<j<n) 
have  joint  density 

. a—.  U-y)w(i-wr-j  _ _ _ 


Pu+uj+iU P,  z,  w)  =  n! 


(i-1)!  0-/-2)!  (n  -  1  —7)!  * 
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We  will  need  their  double  sum 


w-3  **-t 


I  I  Pu+\jj+t(x,  y,  z,  w)  =  n\ 

t- I  y  - 1  +  2 


V3  "y1  l*-y)k  (l-^)"'3"'^ 

*“,(/-!)!  *“0  k\  (n-3-  i-  k)\ 


(2.8) 


i  "y  0- w+z-y)n-}-‘ 

n  ,-i  ( f  —  1 ) !  ( n  —  3  —  i ) ! 


=  n(n  -  l)(n -2)(n -3)(1  -  w  +  z -y  +  x)"  4. 


Proof  of  Theorem  1.  The  expectation  in  E[I(fX)~  I„(fX)]2  is  with  respect  to 
both  the  random  samples  {/JT-i  and  the  random  process  {X(/),0SfSl}  which  are 
mutually  independent.  Performing  first  the  expectation  with  respect  to  the  random 
process  X(t)  we  find,  with  M(t,s)=f(t)R(t,s)f(s), 

E[l(fX)-In(fX)f=E  (IT  M(t,  s)  dt  ds 


-i 

i=  0 


i: 


[M(/,f1)  +  M(/,/i+1)]  dt  -  (ti+l-t,) 


(2.9) 


i  i  tj)  +  M(ti,  tJ+i)  +  M(ti+t,  tj) 

4  j= o y« o 


+  M(fj+l,t/+1)](fl+i-fj)(fy+1-//)J. 

The  interchange  of  expectation  and  integration  is  justified  from 

£{lo l/(0*(,)l  d'}2  =  r  r  l^(0y(^)|JE{|A:(f)A:(s)|}  dtds 


S{1  l/(')l/?l/2(u)d'} 


<00 


since 

£2{|X(t)X(s)|}S  E{X2(t)}E{X2(s)}  =  R(t,  t)R(s,  s). 

The  more  restrictive  condition  Jo/2(f)£(f,  t)  dt  <oo  will  be  needed  for  the  finiteness 
of  the  expected  value  with  respect  to  the  random  samples  and  the  further  interchange 
of  expectation  and  integration.  As  this  is  plain  from  the  following  expressions,  it  will 
not  be  discussed  further. 

To  facilitate  the  computation  of  the  expected  value  with  respect  to  the  random 
sampling  points  /,,*••,  t„,  we  split  the  double  summation  in  (2.9)  into  its  diagonal 
part,  the  part  corresponding  to  the  two  immediate  parallels  to  the  diagonal  and  the 
rest;  which  in  view  of  the  symmetry  of  M(t,  s)  can  be  written  as 

n  n  n  «-l  n-2  n 

11=  I  +2  1  0  =  i+l)  +  2  I  I  . 

1-0 j-0  t—J—O  1-0  i—0  j=  i+2 

Omitting  for  simplicity  the  terms  in  brackets,  which  are  evident  from  (2.9),  we  write 
the  mean-square  error  in  the  form 

£[/(/X)-/b(/X)]2=  f  f  *  M(f,  s)  dtds 
Jo  Jo 


(2.10) 
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(AA.) 


(4A2) 


(A A3 ) 


(AA4) 


-iff  [•••]*  a/, 

i=0  Jo 


+  l-  I  E[--.]A/,A0 

4  i  as  j  =  Q 

Z  ^['  ‘  ’3u=i+nA<iA/i+i 
Z  i  =0 

+  fl  £  £[-*-]A/<Af, 

Z  i=o  j-i+2 


with  A/,  =  /i+)  — fj.  We  now  evaluate  separately  the  cross  term  A,,  the  diagonal  term 
A2,  the  second  diagonal  term  A3,  and  the  off-diagonal  term  A4,  all  clearly  identified 
on  the  left  margin  in  (2.10). 

The  cross  term  A,.  Since  fo  =  0  and  /„+,  =  1,  we  isolate  the  first  and  last  terms  in 
the  sum  £"=0,  and  write 


£  P 

Jo 


[M<\ 0)  +  M(t,  /,)]  dt •  r, 


+  I  £  [M(t,  tj)  +  M(t,  *i+l  )]  dt  ■  (fI+,  -  fj) 


+  £  [M(f,  th)  +  M{t,  1)]  dt-  (1  -f„) 


dt[M(t,  0)  +  M(t,  x)]xp,(x) 


■  [  [  dxdy  [  dt[M(t,x)+M{t,y)](y-x)l  £  pu+i(x,y)  ■ 

J  Jo<x<v<l  Jo  v i  —  1 

f  dj  (  d/[M(f,y)  +  M(f,  1)](1 -jOpnOO- 

Jo  Jo 


We  now  use  the  expressions  in  (2.1)  and  (2.2)  to  write: 


,  =  n  f  dt  J  dx[M(t,0)+M(t,x)]x(\-x)"-' 

Jo  Jo 

+  n(n-l)  I  dt  I  J  dxdy[M(t,x)+M(t,y)](y-x)(l-y  +  x)n 

Jo  J  J0<x<y<\ 

+  n  f  dt  (*  dy[M(t,  y)  +  M(t,  1 )]( 1  -y)yn~> 

Jo  Jo 


dy[M(t,  y)  +  M(t,  l)](l-y)y" 


and  then  we  evaluate  all  inner  integrals  that  can  be  computed  to  obtain 

— A,— — — —  f  M(t,0)dt  +  n  [  j  M(t,x)x(l-x)n~,dtdx 
n  A  1  Jo  Jo  Jo 

+1 1  dldy} 


rt  f  r  M(t, 

Jo  Jo 


y)U-y)yn  ’dtdy+—^—[  M(r,i)dr. 
M  +  1  J0 
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Rearranging  the  terms  we  put  Ax  in  its  final  form: 


(2.11) 


l,  = - — -  f  [M(x,l)  +  M(0,l-x)]dx 

n  +  1  Jo 

-f  I  M(x,y){2-yn-(l-y)n}dxdy. 
Jo  Jo 


The  diagonal  term  A2 .  We  proceed  as  for  /4,.  We  first  separate  the  terms  with 
ro  =  0  and  f„+1  =  1, 

4A2=£{[M(0,0)  +  2M(0,  +  /,)]#?} 

+  X  £{[M(<l,f1)  +  2M(/„/I+l)+M(r,+1,f1+1)](f,+1-r,)2} 

<*1 

+  £{[M(t„,  i„)  +  2M((„  1)  +  M(l,  l)](l-/„)2} 
and  then  using  (2.1)  and  (2.3)  we  obtain 


2  =  ”  J" 

Jo 


4A2=  n  [M(0)0)  +  2M(0,x)+M(x,x)]x2(l-x)"'1  dx 


0<x<y<i 


[M(x,  x)  +  2M(x,  y)+  M(y,  y)] 

11 

•  (y-x)2(l->'  +  x)"~2  dxdy 


-n  f  [M(y, 
Jo 


y)  +  2M(y,  1)  +  M(1, 1)](1  -y)2yn~'  dy. 


We  now  evaluate  all  inner  integrals  that  can  be  computed  and  regroup  terms  to  reach 
the  final  expression 

A^^rf [,[m(x,i)+m(o)i-x)](i-x)v-1* 

2(n  +  l)(n  +2)  2  J0 

(2.12)  +  '  M(u,  «)  [— ^  +  rr— [Mn+I  +  (1  ~  M)"+1]  [“"  +  (1  -  “)”]}  du 

Jo  in  +  1  2(n  +  l)  2  J 

+  \n(n-\)  ||  M(x,y)(y-x)2(\-y  +  x)n-2dxdy. 

J  Jo<x<y<r\ 

The  second  diagonal  term  A}.  We  first  split  off  the  sum  X  "Jo*  the  (>rst  and  last  terms 
involving  to  =  0  and  tn+,  =  1, 

2 Aj  «  E{[M(0,  t,)  +  M(0,  t2)  +  M(t, ,  t,)  +  M(t, ,  t2)]t,(t2- /,)} 

n—2 

"f"  X  <(+i)  +  M(tl%  f,+2)  +  M(ti+x ,  fj+i)  +  M (tj+i ,  fj+2)] 

l  *  I 

'  (*f  +  l  —  0(* i+J  —  f(+l)} 


TRAPEZOIDAL  MONTE  C  ARLO  I NTECj RATIO N 


237 


Next  we  use  the  values  of  the  bivariate  densities  p,.;  and  p„  ,  „  from  (2.3)  and  the 
sum  of  the  consecutive  bivariate  densities  in  (2.5)  and  write 


2/4,  =  n(n 


[M(0,  x)  +  M(0,  >’)  +  M(x,  x)  +  M(x,  y)] 


-„}} 

J  J  0-  X  *  V  ■ 

•  x(y-x)(l  - y)n ~2  dxdy 

+  n(n-  l)(n  -2)  j"J|  [M(x,  y)+  M(x,  z)+  M(y,  y)+  M(y,  z)] 

•  (y  —  x)(z ->’)(!  -z  +  x)"  '  dxdydz 


+  n(n-l) 


[M(x,y)  +  M(x,  \)+M(y,  y)+M(y,  1)] 
•  (,v-x)(l -y)xn~2  dxdy. 


We  again  evaluate  all  inner  integrals  which  can  be  computed  and  regroup  terms  to 
derive  the  following  expression,  after  considerable  algebra, 


(2.13) 


2(n  + 1)  Jo 


[M(  x,  1)  +  M(0, 1  —  x)]xn_2(l  —  x)  | 

M(u,  u)[l-un+,-(l-u)"+']du 

M(x,y)(y-x){--^[x"-'  +  (l ->>)"'] 


dx 


0 <x<><] 

+  ^(n-2)(>’-x)2(l  —  y  +  x)n_3H  — —  (1— ^  +  x)"_1|  dxdy. 
o  n  —  l  J 

The  off-diagonal  term  A4.  We  first  isolate  the  terms  involving  the  points  fo  =  0  and 
f„+i  =  1  by  splitting  the  double  sum  into 

n-2  n  n - 1  n-2  n-3  n - 1 

I  I  =(i  =  0,j=n)+  l  (i  =  0)  +  l  (j  =  n)+  £  £  . 

i=0  j=i+2  7  =  2  i=l  i  =  l  7=1  +  2 

We  thus  write  A4  as 

2/44=E{[M(0,  /„)+M(0,l)+M(r1,/n)  +  M(r1,i)3rI(i-rn)} 

+  "i  E{[M(0,  (;)  +  M(0,  o+1)  +  M(/„  tj)+  M{tu  />+1)]f,(o+,-  tj)} 

7-2 

+  £  £{[M(f„  /„)  +  M{tit  1 )  +  M(ti+1 ,  r„)  +  M(/,+1, 1 )](/,+,-  tf)(l  -  f„)} 

i-1 

+  "l  "l  £{[M(/f,o)  +  M(t(,o+I) 

/=*  1  >»i+2 


+  M(tl+ 1 ,  (7)  +  ,  fj+i)](f(+,  —  h)(tj+i  ~  tj)}- 


1 
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In  calculating  the  expectation  we  u»e  Tor  the  first  term  the  value  of  pl  m  from  (2.2),  for 
the  second  term  the  sum  in  (2.6),  for  the  third  term  the  sum  in  (2.7),  and  for  the  fourth 
term  the  sum  in  (2.8),  and  obtain 

2A4=n(n-\)  f  [  [M(0, y)  +  M(0, !)+  MU y)  +  M(x,  I)] 

•  x(l -y){y-x)a~2  dxdy 

+  n(n  -  1 )( ft  -2)  I  I  I  [M(0, y)  +  M(0,  z)  +  M (x,  y)+  M(x,  z)] 

•x(z-x)(l-z  +  y-x)"  'dxdydz 

+  n(n  -  l)(n  -2)  f  f  f  [M(x,z)  +  M(x,  1)  +  M(y.  z)  +  M(y,  1)] 

J  J  Jo*  *•  yv:*  | 

•  (y  -  x)(l  -  z)(z  -y  +  x)"  'dxdydz 

+  n(n-l)(n-2)(n-3)  JJJJ  [M(x,  z)+  Mix,  w)  +  M(y,  z) 

+  M(y,  w)](y  -  x)(  w  -  z)(l  -  w  +  z  -y +  x)"  4  dxdydz  dw. 

Now  we  evaluate  all  inner  integrals  that  can  be  computed  and  we  regroup  similar 
terms.  After  extensive  but  routine  calculations  we  find 


(2.14a) 


*4=  [  f  Mix, 

Jo  JO 

•{— - [x"  +  (l-x)"]+—  [x"+,  +  (l-x)"+,]  +  x(l-x)" 

In  + 1  n  +  1 

1  ,  w  J  xn~2  x"-1  x"  X-  11  , 

2  L  3(rt  -2)  n  - 1  n  3(n+l)JJ 

+in(n-l)[[  M(x,y)A(x,y)  dxdy 

£  J  Jo<x<y< 1 


where 


A(x,  y )  =  x(  1  -  y  )(y  -  x)-2+ x  {2  ( I  -  x)[(y  -  x)" 


+(i-yr-2]+^-7[(y-xr-’+(i-yr 

fl  —  1 


+  0 


-4£- 


->} 


y[(>'-x)n~2  +  x’’-2]  +  ~ — 7  [(y  -x)"'1  +  x” 
n  —  1 


■>} 


(2.14b)  +— [y"  +  (l  -x)"  +  (l  —  y  +  x)"]— —  {x"  +  (l  -y)"  +  (y  — x)"] 
n  n 


n  —  l 
n-2 
n-1 


[y"",  +  (l-x)-,  +  (l-y  +  x)"-'] 
[x"-'  +  (l-y)"-,  +  (y-x)"-'] 


+  x"~2y(l  —  y  +  x)  +  (l  -y)"”2(l  -x)(l  -y +  x)  +  (y  -x)"_2(l  -x)y 


+  J<«- 2)(» 


-3>{ 


(1-y  +  x)"  (1-y  +  x)"  1  (1  — y +  x)"~2  (1  - y  +  x ) 


3n 


n-1 


n-2 


3(n 


till!] 

-3)  J 


?;-V 

A IL  \  . 


-4 
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We  now  substitute  in  (2.10)  the  expressions  for  <4,  to  A4  we  derived  in  (2.11)  to 
(2.14),  and  after  grouping  similar  terms  and  some  algebra  we  arrive  at  the  expression 
in  (1.6).  □ 

3.  The  rate  of  convergence.  In  this  section  we  determine  the  rate  of  convergence 
to  zero  of  the  mean-square  error  given  in  Theorem  2. 

We  will  use  the  following  expression  for  the  integral  for  a  function  F  with  i.  ( =  1 ) 
continuous  derivatives: 

f  L  <— n*'1  (-1!1  P 

(3.1)  F(x)x"<fa=[4rF,M,(l)  +  LlF  Fu\x)x"'ldx 

Jo  /-I  W  n  Jo 

where  n  is  a  positive  integer  and 

(3.2)  n,,l  =  (n  +  l)  •••  («  +  /). 

Expression  (3.1)  is  obtained  by  repeated  integration  by  parts  or  use  of  Taylor's 
expansion  for  F(x)  about  the  point  1. 

We  will  also  use  the  following  version  of  the  approximation  of  a  function  by  a 
delta  sequence  (whose  proof  is  standard). 

Lemma.  Let  F(x)  and  {K„(x)}*.t  be  Borel  functions  defined  on  [0, 1J.  If  F  is 
bounded  and  (left)  continuous  at  1  and  if  the  kernels  K„  satisfy  the  following  conditions: 


(i)  f  Kn(x)dx=  1  for  all  n, 

Jo 

(ii)  f  |K„(x)|  dxS  C<oo  for  all  n, 

Jo 

(iii)  lim  [  \K„(x)\dx  =  0  for  all  8  e  (0,  1), 

Jo 


lim  F(x)Kn(x)dx  =  F(  1-). 


In  particular,  when  F  is  as  in  the  lemma  we  have 


(3.4a) 


(3.4b) 


(3.4c) 


lim(n  +  l)[  F(x)xndx  =  F(  1-), 

J0 

lim  (n  +  l)(n  +  2)  [  F(x)x"(l -x)  dx  =  F(l-), 

"-00  Jo 

lim  ^;(n  +  l)(n  +  2)(n  +  3)  |  F(x)x"(l  - x)2  dx  =  F(l-). 
"-°°2  Jo 


Proof  of  Theorem  2.  We  proceed  by  determining  first  more  detailed  expressions 
in  inverse  powers  of  n  for  each  of  the  integral  terms  in  (1.6),  which  are  then  combined 
to  produce  the  rate  of  convergence  of  the  mean-square  error.  We  denote  by  B ,  the 
first  integral  (1.6b),  by  B2  the  second  integral  (1.6c)  and  by  B,  the  third  double  integral 
(1.6d).  For  convenience  we  set  M(x,  y)  =f(x)R(x,  y)f(y). 
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The  sectional  integral  B,.  Putting  F(x)  =  Af(x,  1)  +  M(0, 1-x)  and  using  (3.1) 
with  L  =  2  (and  changing  variables  to  y  =  1  -  x  for  the  final  term)  we  obtain  after  some 
algebra 


F(0)  F(l)  F'(0)  F'(l) 
°'~~2  na)  n(2)  2nw  2n(J> 


F"(l  -x)x"+3  dx 


+  3n  -2)x 


n+3 


l)(n  +  2)(»  +  3)  2(n  +  2)  4( 


n  +  l)J 


It  is  easily  checked  that  the  polynomials  kn(x)  within  braces  in  the  last  integral  are 
positive  on  (0, 1 ),  with 

1  Ux)dx  =  ^’ 

and  that  the  kernels  Kn(x)  =  2n<4)M*)  satisfy  the  conditions  of  the  lemma.  It  then 
follows  from  (3.3)  and  (3.4a)  that 


F(0)  +  2F(1)  F'(1)-F'(0)  F"(l)  -  F"(0)  ,  ,  _4% 
- 2*“"  *  2n">  +°<n  >■ 


and  using  the  form  of  F  to  express  the  coefficients  in  terms  of  M  we  find 
B,  =  -^{M(0,0)  +  M(ltl)  +  M(°,l)} 


(3.5) 


+^{M'-0(1,  l)-M,  o(0,0)-M*  °(0,  l)  +  Mo-'(0, 1)} 


2n 


(4) 


{M2O(0, 0)  +  Af20(l,  1)  -  M  2O(0, 1)  -  Afo,2(0, 1)}  +  o(n-*). 


The  diagonal  integral  B2 .  The  integral  B2  may  be  written  in  the  form 

B2  —  , ,  3,  tt  [  M(u,u)du+\  {M(u,  u)  +  M(l-u,  1  —  m)} 

2(n  + 1)  Jo  Jo 

l2(n  +  l)  2  J 

Since  F(u)  =  M(u,  u )  +  M(  1  -  u,  1  -  u)  has  two  continuous  derivatives  we  obtain  from 
(3.1)  with  L  =  2,  applied  to  the  second  term. 


M(u,  u)  du 

3 

2n<2) 

un+2l 

1 2n<3) 

~2n<2,j 

'(1) 


The  polynomials  Mu)  within  braces  in  the  last  integral  are  negative  on  (0, 1)  with 

5 


[' 


M«)  du  =  - 


2n(4” 
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and  it  is  easily  checked  that  the  kernels  K„  =  (2 /5)n<4)k„  satisfy  the  conditions  of  the 
lemma.  It  then  follows  from  (3.3)  that 


=  [  M(u,u)du- 

2(n  +  \)  Jo 


3F(1)  2F'(1)  5F"(1)  ,  ,  _4x 
- oT - °  n  ' 


2n 


(2) 


n 


and  using  the  form  of  F  we  obtain 
3 


B?  = 


2(n  + 1) 


f  M(u,u 

Jo 


2  n 


(3.6) 


+4t{W-°(U1)-m,-°(0,0)} 

n 


.<*) 


{M2’°{  0, 0)  +  M2-°d,  1)  +  Mw(0,0)  +  Mm(1,1)}  +  o("_4). 


The  double  integral  B3.  In  the  expression  of  B3  in  (1.6d),  by  changing  variables 
appropriately  in  each  of  the  six  double  integral  terms  and  grouping  separately  the  first 
three  and  the  last  three  terms,  we  can  write  B3  in  the  following  form: 


B3 

(3.7)  + 


•nil' 

ill; 

=  [  F(x)x"dx+[  G(x)g„(x)  dx 
Jo  Jo 


[3M(x,  y)  +  3M(l— x,  1— y)  + A/(y  — x,  y)]  dy 

1 


}*" 


dx 


M(x-y,l-y)  dy 


}|i(n-2)(M 


+  3)xn  — -  n  x  + 


~n(n-l)x"-2J 


dx 


with  the  obvious  identification  for  F,  G,  and  g„.  Since  /  is  assumed  to  have  only  two 
continuous  derivatives,  so  do  F  and  G.  We  therefore  use  first  (3.1)  with  L  =  2  and 
then  identify  those  terms  in  F",  G"  which  are  not  differentiable  for  separate  treatment. 
We  first  obtain  (after  some  algebra) 


„  F(l)  F'(l)  .  1  f  ,  „+2  , 

B  — - - F  (x)x  dx 

n  +  1  n<2>  n<2)J0 

(3-8)  .....  .....  .. 


3G(1)  +G 


2(n  +  l) 


m+  r  <riX)ii*z _ «_x~.+ix.i 

n<2)  Jo  °Wl  4n<2)  X  2(n  + 1) X  Vi** 


Evaluating  2F(1)-3G(1)  and  G'(l)-F'(l),  and  denoting  by  h„(x)  the  polynomial  in 
braces  in  the  last  integral,  we  have 


(3.9) 


P  F"(x)xn+2dx+  f  G"{x)h„(x)  dx. 
n  Jo  Jo 


From  the  definition  of  F  and  G  in  (3.7)  we  evaluate  their  second  derivatives  after 
which  we  separate  those  terms  which  are  not  differentiable.  We  thus  write 


F"  =  F,  +  F2>  G"=G,  +  G2 


where 


(3.10a)  F,(x)  =  -1  [M'-°(x,  x)  - M,  0(  1  -x,  l-x)]-&M°'(0,  x) -  M,  0( 0,  x)], 
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(3.10b)  F2(x)  =  ^j  {3M20(x,y)  +  3M2'°(l-x,l-y)+M20(y-x,y)}  dy, 

(3.1  la)  G,(x)  =  -M°'( 0, 1  - x)  +  M'-0( 0, 1  - x), 

(3.11b)  G2(x)=  I  M2  0(x-y,\-y)  dy, 

Jo 

and  F,,  G,  are  differentiable,  while  F2,  G2  are  not.  Proceeding  as  before,  using  (3.1) 
with  L  =  1,  we  obtain 


B 


f  F,(x)x"+2  dx+  f  G,(x)h„(x)  dx 
n  Jo  Jo 


(3.12)  =^--4f,F'I(x)x"+3dx 

n  n  Jo 


+  G(l)x{0} 


~2)(n  +  3)  n+i 


4  n<3>  *  2n(2> "  4(n- 

It  is  easily  checked  that  the  polynomials  kn(x)  within  braces  in  the  last  integral  are 
positive  on  (0, 1 )  with 

3 


;x"+2  +  - 


— 7X"  +  ,]dX- 

1  +  D  J 


r 


k„(x)  dx  = 


2n<4) 


and  that  the  kernels  Kn(x)  =  (2/3)w<4,kn(x)  satisfy  the  assumptions  of  the  Lemma. 
Thus  from  (3.3)  and  (3.4a)  we  have 

f,(d  f;u)  3g;q) 


b3i  = 


.(3) 


Evaluating  F,(l)  and  2Fj(l)  +  3G',(l)  from  (3.10a)  and  (3.11a)  we  finally  find 
®j..  =r~(5j{9(M,o(0, 0)-M1-O(l,  l)]  +  Ml,o(0,  l)-Mo,(0, 1)} 

2ft 


(3.13)  +^(9M2'°(1, 1 )  +  6M2,o(0,  0)  +  Mo,2(0,  1) 

+9Ml,l(l,  1)  +  12M,,I(0, 0)-M'  '(0,  l)}  +  o(n-4). 

To  complete  the  evaluation  of  B ,  in  (3.9)  we  need  to  evaluate 

i?3.2  — — 72)  f  F2(x)xn+2dx+  f  G2(x)hn(x)  dx. 
n  Jo  Jo 

While  F2  and  G2  are  not  differentiable,  when  M(t,  s)  =f(t)R{t,  s)f(s)  is  substituted 
in  (3.10b),  (3.11b),  some  of  the  terms  in  the  resulting  expressions  are  differentiable 
and  we  first  deal  with  these.  We  thus  decompose  F2  and  G2  into 

f2  =  f2i, + f2-2  ,  g2  =  g2>,  +  g2i2 

where 

F2.,(x)  =  X-  J'  {6 f(x)R(x,  y)f{y)  +  6f(l  -x)R{l  -x,  1  -y)f(  1  -y) 

(3.14a)  +  2f(y  -  x)R  ifi(y  -  x,  y)f(y)  +  3 f(x)R2-°(x,  y)f(y) 

+  3/(l-x)R20(l-x,l-y)/(l-y) 
+f(y-x)R2-°(y-x,y)f(y)}dy, 
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(3.14b)  F2.2(x)  =  i  J‘  {3r(x)R(x,y)f(y)  +  3r(l-x)R(l-x,  1  -y)f(  1  -y) 

+/*(>’  -x)R(y~x,  y)f(y))  dy, 

(3.15a)  G2Jx)  =  £  {2f(x-y)R'  °(x-y,  l~y)f(l-y) 

+f(x-y)R20(x-y,l-y)f(\-y)}  dy, 

(3.15b)  G2.2(x)  =  [X  nx-y)R(x-y,\-y)f(\-y)dy. 

Jo 

As  F2 1  and  G2,i  are  differentiable,  the  term 

#3.2.1  -  f  F2't(x)x"+2  dx  +  [  G2A(x)hn(x)  dx 
n  Jo  Jo 

can  be  evaluated  by  using  (3.1)  with  L=  1  to  obtain 

*3A,-“ Bi^iO)— ^  f  F2  ,(x)x"+1  dx  + G2  ,(1)  x  {0}-  [  Gi.,(x)k„(x)  dx 
n  n  Jo  Jo 

with  the  same  polynomials  k„(x)  as  in  (3.12).  It  then  follows  as  for  B, ,  that 


B 


3.2.: 


_F2,,(1)  Fiji)  3G2jI(1~) 


„<3) 


«<«) 


+  o(n'4) 


n"  n'”  2n 

and  evaluating  F2  ,(1)(=0)  and  2F2  ,(1)+3G2>I(1-)  from  (3.14a)  and  (3.15a)  we  obtain 
B3.2.,  =  ^{6B,  0(1,  1)/(1)/'(1)  +  2R,  o(0,  l)/'(0)/(l) 

+  3R20(1,  1)/2(1)  +  R20(0,  l)/(0)/(l) 

U'l6;  -3  j1  [2B10(u,  u)f(u)f"(u)  +  3R2’°(u,  u)f(u)f(u) 


+  /?3,0(«-,  u)/2(u)]d«J  +  o(n-4). 

To  complete  the  evaluation  of  B3  2  we  finally  need  to  evaluate 

#3,2.2 A -75)  f  F2>2(x)xn+2dx+ I  G22(x)/i„(x)  dx 
n  Jo  Jo 

where  F2i2  and  G2  2  are  given  in  (3.14b)  and  (3.15b)  and  h„  is  the  polynomial  in  braces 
in  (3.8).  Substituting  F2 ,2  and  G2>2  and  isolating  the  nondiiferentiabie  factor /"  we  can 
write  it  in  the  form: 

#3.2.2  «  | '  /"(x)x"+2 1 1 ' '  R(x,  y)f(y)  dy }  dx 

+ T  r(  1 "  x)x"+2  { I R  ( 1  -  * 1  ■  * )/( 1  ■  y >  dy}  dx 

*(u,o)/T»)(o-«)"+2di,j  rfu 
+  |  /"(«){J  #(u,t>)/(t')Ml  +  w-t>)dfj  dw. 
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We  denote  the  four  terms  by  7],  i  =  l,  2,  3,  4,  and  the  corresponding  functions  in 
braces  by  Ht.  For  the  first  two  terms  we  use  the  Taylor  expansion  about  1, 

H(  x)  =  H(  1 )  ■ -  H'(  1 )( 1  -  x)  +  i  H"(  tx )( 1  -  x)2 

where  the  intermediate  point  tx  belongs  to  (x,  1)  and  depends  continuously  on  x.  We 
thus  find 


T,  +  Tz~  2n<2) 


(*)*"+2(l  ~x)dx+-  f”(x)H”l(tl'X)xn+2(l-xy  dx 


{  J '  /"(*)*"+2(  1  -  x)  dx  + X-  £  f'( 

-HUD  f  f"(l-x)xn+2(l-x)dx 
Jo 

+\  JV"(l-x)//J(r2.Jx"+2(i-x)2  dx] 


and  from  (3.4b),  (3.4c)  we  obtain 


r, + r2 {-//;(1)/"(1)  -  H  i(D/”(0)}+ o(«‘4)+ o(n-5) 


and  finally 


(3.17)  r,  +  r2  =  —fj;  {R(  1 , 1 )/( 1 )/"( 1 )  +  R(0, 0)/(0)/"(0)}  +  o(  n ~4). 

In 

For  the  third  term,  T3,  we  first  integrate  by  parts  the  expression  of  Hy  to  write  it  in 
the  form 

H3(«)  =  4-/?(M)/(1)(l-«rJ - “T  {'  Dv[R(u,v)f(vmv-u)n+idv 

n  + j  n  + i  Ju 

where  D„  denotes  partial  derivative  with  respect  to  v,  and  substitute  in  T3 : 
r3  =  |o'  /"( u ) « ( «, 1 )  ( 1  -  « ) " +3  du 

-A)  f  f  f"(u)Du[R{u,  v)f{v)](v-uy+}  dv. 

J  Jo<u<v<\ 

Since  the  integrand  of  the  double  integral  is  bounded  and  JJu<l)(o-u)kdudt;  = 
[(fc  +  l)(fc  +  2)]_1,  the  double  integral  is  0(n-5)  and  thus  by  (3.4a)  we  obtain 


(3.18) 


T^^rjwrWR^D+oin-4). 


For  the  fourth  term  T4  we  likewise  first  integrate  by  parts  the  inner  integrand  HA. 
Noting  that  h„  =  k’„  where  k„  is  the  polynomial  within  braces  in  the  integral  in  (3.12), 
we  obtain  by  integrating  by  parts 


/f4(u)=  [  R(u,  l  +  u-z)f(l  +  u-z)k'„(z)  dz 
J  U 

=  /?(«,  u)f(u)k„(l)-R(u,  1)/(1  )k„(u) 


-l'o, 


[R{u,l  +  u-z)f(\  +  u- z)]k„{z)  dz. 
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Since  k„(  1)  =  0,  and  as  was  pointed  out  following  (3.12),  (2/3 )n'*'kn  satisfies  the 
assumptions  of  the  Lemma,  we  obtain  by  (3.3), 


(3.19) 

r. 


(z)R(z, !)/(!)  + 


(J/c„(z)  dz 


3 

2/i<4) 


|  f"(u)D2[R(u,  l  +  u-z)/(l  +  u-z)]  duj 

|/"(1)H(1, 1 )/( 1 )  -  JJ /"( w)[ /20-1  ( m,  u)f(u)+R(u,  u )/'(«)]  duj  +  o(n-4). 
Putting  together  the  expressions  in  (3.17)  to  (3.19),  we  find 
B3, 2.2  =  ^  { : 3  ■ *  (°.  0)/(0)/"(0)  +  R  (0, 1  )/"(0)/(  1 ) 

+  3  J  f"(u)[R°\u,  u)f(u)  +  R(u,  u)f(u)]  du}  +  o(n~4). 

Now  from  (3.9),  (3.13),  (3.16),  and  (3.20),  we  derive  the  final  expression  for  B, : 
B}  =  ~2(7m)  1  M(“’  U)  du+~n^  {2M(°.0)  +  2M(1, 1)  +  ^  M( 0, 1)} 


(3.20) 


+  ^{9[M'-O(0,0)-M1O(l,  1)]  +  M'  °(0, 1)-  Mo,l(0, 1)} 


(3.21) 


+  ^jJ6M2O(0,  0)  +  9M2-0(l,  1)  +  Mo,2(0,  1)  +  12Mij(0,  0) 

+  9M1-1(1, 1)  — Af'^O,  1)  +  3R(0, 0)/(0)/"(0)  +  /?(0,  l)/"(0)/(l) 
+  6B,  0(1,  l)/(l)/'(l)  +  2B,  o(0,  l)/'(0)/(l)  +  3/?20(l,  1)/2(1) 

+  R20( 0,  l)/(0)/(l)  +  3  I  [/''B/'-/"R,  0/-3/'/?20/](M,  u)  du 

Jo 

-3  £  K3->-  u)f\u)  du}  +  o(n-4). 


The  rate  of  convergence.  Substituting  the  expressions  (3.5),  (3.6),  (3.21)  of  B},  B2, 
B3  into  the  expression  (1.6)  of  the  mean-square  error,  we  find  that  all  lower-order 
terms  cancel  and  we  obtain  (1.7),  where  the  constant  C(f  R)  is  readily  identified  from 
the  coefficients  of  (n<4))-1  in  (3.5),  (3.6),  and  (3.21).  In  addition  to  expressing  M  in 
terms  of  /  and  R  in  the  coefficient  of  (n(4))~'  in  (3.21),  we  also  use  the  following 
expressions  for  some  of  the  integrals  involved,  which  follow  by  integration  by  parts: 

jo  [/"«/'](«,  «)  du  =  ± {[/'R/'](l,  l)-[/'K/'](0,0)}- J’  [/'Rl0/'](u,  u)  du, 

f  [fnRi  0f](u,  u)  du  =  [/'«' °/](l ,  l)-[/'B,  o/](0,0) 

Jo 

-  P  [f(R2'°f+ R'-'f+ R'-°mu,  U)  du. 

Jo 

The  resulting  expression  of  the  asymptotic  constant  C(f  R)  is  given  in  (1.8).  This 
completes  the  proof  of  Theorem  2.  □ 
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The  asymptotic  constant.  The  expression  of  C{fR)  given  in  (1.8a)  cannot  be 
symmetrized  any  further  under  the  current  assumptions  (note  the  lack  of  symmetry  in 
the  constant  terms  involving  /?'•'  and  in  the  integrals).  However,  if  R(t,  s)  is  further 
assumed  to  have  continuous  mixed  partial  derivatives  of  order  3  throughout  the  unit 
square  (rather  than  olf  its  diagonal — as  has  so  far  been  assumed),  then  by  integrating 
by  parts  we  find  (using  the  obvious  shorthand) 

(3.22)  f  (Ru-2R20)jr  =  12[(R,,,-2R2-0)f2]o+  f  R30/2 

Jo  Jo 

and  thus  the  integral  terms  in  C(f  R)  can  be  evaluated.  This  produces  the  following 
symmetric  expression: 

2 Csym(f  R)  =  iR(0,0)f’(0)2+iW,  1)/'(1)2-K(0,  l)/'(0)/'(l) 

+  R,o(0,0)mf(0)  +  R'-°(l,  1)/(1)/'(D 

(3-23) 

~RU0( 0,  l)/(0)/'(l)-K°-‘(0,  l)/'(0)/(l) 

+  j/?l’1(0, 0)/2(0)  +5R1,i(1,  1)/2(1) -  /?m(0,  l)/(0)/(l). 

It  is  easily  checked  that  this  can  be  written  as  in  (1.10),  so  that  under  these  slightly 
more  stringent  assumptions,  which  fall  short  of  guaranteeing  two  quadratic-mean 
derivatives  for  X,  we  have  (1.11).  In  view  of  (3.22)  and  (3.23),  the  general  form  of  the 
asymptotic  constant  in  (1.8a)  can  be  written  as  in  (1.8b). 
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